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Abstract 

We study the role of unsteady lift in the context of flapping wings in birds flight. Both aerodynamicists 
and biologists attempt to address this subject, yet it seems that the contribution of the unsteady lift still 
holds many open questions. The current study deals with the estimation of unsteady aerodynamic forces 
on a freely flying bird through analysis of wingbeat kinematics and near wake flow measurements using 
time resolved particle image velocimetry. The aerodynamic forces are obtained through unsteady thin 
airfoil theory and lift calculation using the momentum equation for viscous flows. The unsteady lift is 
comprised of circulatory and non-circulatory components. Both are presented over wingbeat cycles. Using 
long sampling data, several wingbeat cycles have been analyzed in order to cover the downstroke and 
upstroke phases. It appears that the lift varies over the wingbeat cycle emphasizing its contribution to 
the total lift and its role in power estimations. It is suggested that the circulatory lift component cannot 
assumed to be negligible and should be considered when estimating lift or power of birds in flapping 
motion. 


Introduction 

Interest in the unsteady aerodynamics of flapping wing flight has been rekindled with the development 
of micro air vehicles (or MAVs). These MAVs fly at low Reynolds numbers, where many complex flow 
phenomena take place within the boundary layer. For example, separation, transition, and reattachment 
(of the boundary layer) can occur within a short distance along the surface of the wing and can dramatically 
affect the performance of the lifting surface. Despite these challenges engineers are not without prior 
information because nature has produced numerous examples of biological flying machines that have 
evolved over millions of years to efficiently fly in the low-Reynolds-number regime. One such example 
is the flapping flight mechanism in which the wings are not only moving forward relative to the air, 
but also flap up and down, bend, twist and sweep, resulting in a complicated unsteady-aerodynamic 
motion. Understanding the role of unsteady flapping flight will help in designing more efficient micro-flying 
vehicles [?]. 

The flapping motion associated with the unsteady effects generally leads to enhancement of bound 
vortices on the lifting surface, which eventually detach, convect into the wake, and interact with other 
vortices [?]. Due to the interaction between the bound vortices on the lifting surface and the vortices in 
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the wake, the performance of an unsteady wing is coupled with the formation and distribution of vorticity 
shed throughout the wings cycle of oscillation [?,?]. 

A useful theory to approximate unsteady aerodynamic loads is the unsteady thin airfoil theory. The 
roots of the theory were originally developed by Glauret [?], who considered simple harmonic motion. 
However, the complete solution of estimating the time dependent loads for a harmonically oscillating 
airfoil in inviscid, incompressible flow was first published by Theodorsen [?]. Tlieodorsen’s work was 
further complemented by von-Karman and Sears [?] who simplified the equations and presented the 
general unsteady thin airfoil theory. In addition to simplifying the equations, von-Karman and Sears 
considered also the problem of a thin airfoil moving through a vertical gust field. The effect of unsteady 
inflow conditions on aerodynamic forces is considered in many applications, for example in helicopter 
aerodynamics [?,?]. 

The role of unsteady forces is significant in estimating the aerodynamic performance of birds flight in 
flapping motion [?]. As living organisms, birds are subject to selective pressures, as such, one may assume 
that they operate their wings in a highly efficient manner. This notion is supported by the tendency of 
birds, as well as many other animals, to operate in a limited Strouhal number range between 0.2 and 
0.4 [?,?]. There are many factors differentiating flapping of a bird’s wings from heaving or pitching of two 
dimensional airfoils [?]. These differences include the presence of a body, three dimensionality of the wing 
and its unique motion. Ben-Gida et al. [?] compared the formation of steady to unsteady drag at the near 
wake of a European Starling (Sturnus vulgaris). It was demonstrated that the unsteady drag component 
at the transition stages of the wingbeat phase reduce the total drag. 

To model the time dependent aerodynamic forces acting on a section of a wing it is natural to start 
with a quasi-steady approach. The estimation of lift from the PIV measurements behind flying birds is 
done by applying the classical Kutta-Joukowski theorem, L = pUT, where p is the fluid density, U is the 
wind speed, and T is the circulation calculated from the vorticity fields. Whether the lift is estimated 
from the Trefftz plane or from the streamwise-normal plane, it is assumed to be quasi-steady [?,?]. In 
order to estimate the quasi-steady lift, it is sufficient to capture a portion or an entire wingbeat cycle. 
Former work has shown that circulation can be estimated from a single instantaneous vector map [?] or 
from synchronized velocity maps triggered to match various phases within the wingbeat cycle [?]. Or, 
using several consecutive velocity maps of which a full wingbeat cycle has been reconstructed and lift 
was estimated for a series of velocity fields capturing the far wake behind a freely flying bird [?,?]. Yet, 
incorporating the unsteady effects in lift estimations is lacking. One of the challenges in estimating the 
evolution of lift over time is the need to measure the wake using a technique that introduces high spatial 
and temporal resolution over a large period of time. In the present work, the near wake of a freely flying 
European starling ( Sturnus Vulgaris) has been selected as a case study of unsteady wing aerodynamics [?]. 

The aim of the present work is to evaluate the unsteady sectional lift of a flapping wing, based 
on experimental data acquired in the near wake of freely flying European starling ( Sturnus vulgaris) 
using long-duration high-speed Particle Image Velocimetry (hereafter PIV) [?]. In the case of a flapping 
wing [?], the boundary layer over the wing often experiences an early transition to turbulence due to 
the unsteady motion and remains attached for higher angles of attack (compared to airfoil in steady 
condition). Such re-attachment of the boundary layer allows the use of the unsteady thin airfoil theory 
for lift estimation. As a first approximation the wing is estimated as rigid plate undergoing translational 
motion using kinematic relations [?] based on von-Karman and Sears [?] model. Then, using high-speed 
PIV data, the unsteady portion of the lift is estimated over few wingbeat cycles and compared to the 
rigid plate lift approximation. 
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Materials and Methods 


Theoretical modelling 

The unsteady motion of a viscous fluid about a lifting surface is always accompanied with shedding 
of vortices into the wake. The addition point of the unsteady thin airfoil theory derivation is through 
the fact that the circulation around the surface is varying continuously and two-dimensional vorticity is 
shed off into the wake, thus allowing the use of planar wake assumption and ignoring the effect of wake 
deformation. In the far wake, the vortices roll-up under their own self-induced velocities to form complex 
wake patterns. Despite this difficulty, the contribution of the wake vortices is not significant since the 
influence of the shed vortices reduces with increasing distance. However, in the near wake, where the PIV 
measurements were conducted, the flow structures are simplified, thus allowing the use of planar wake 
assumption, ignoring the effects of wake deformation. Here we derive the unsteady thin airfoil theory. For 
clarity, we use the unsteady Bernoulli equation to evaluate the time-dependent loads. This derivation, to 
the best of our knowledge, has not been reported elsewhere. 

We follow the classical definition of the unsteady thin airfoil theory [?,?] notation where the chord is 
equal to c = 2 (see Figure [I]). The leading-edge of the lifting surface is placed at x = — 1 and trailing-edge 
at x = 1, whereas the mid chord is placed at the origin of the coordinate system at x = 0. The y-axis is 
perpendicular to the flow and the z-axis is in the spanwise direction of the wing. 

The analysis is initiated by assuming that the vorticity distribution bound to the wing section j(x, t ) is 
the sum of the quasi-steady vorticity distribution 7 0 (x,f) that would have been produced in quasi-steady 
motion and a wake-induced vorticity distribution 71 (x,t). 

la{x,t) = lo{x,t) + 7 i(ir,t) ( 1 ) 


The total circulation about the airfoil due to both the quasi-steady vorticity distribution and that induced 
by the wake is 

r a = r 0 + Tj (2) 


where 


and 


T 0= J 7o(x,t)dx 
Fi =J 71 (x,t) dx 


( 3 ) 

( 4 ) 


One of the fundamental assumptions in inviscid aerodynamics is that, according to Kelvin circulation 
theorem, the total circulation of a system is equal to zero. The total circulation about the wing section is 
a sum of the quasi-steady circulation To, that would be produced if the total circulation would not have 
been affected by the presence of the wake, and the wake-induced circulation Fi. When the circulation 
around the wing section T a is balanced with the circulation produced by the wake at every time step t 
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J 7o0r,t)dx + 
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dt 


J 71 (x,t)dx + 
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dt 


7 w(€,t) d£ = 0. 


( 5 ) 


Wake-induced circulation 


The effect of the wake vortices is evaluated in accordance with the thin airfoil theory, where the Joukowski’s 
conformal transformation is used to transform a circle in the z' plane to an airfoil at the 0 plane. The 
transformation relating the two planes is 
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A single vortex element defined as T at a distance of £ in the 0 -plane is located in the z'-plane at a 
distance of rj. To create a unit circle which transforms into a flat plate in the z plane, another vortex 
element with an opposite sign has to be introduced inside the unit circle at a symmetric point, which is 
l/? 7 . Thus, along a unit circle the induced velocity magnitude is 


vg 


r' r' 

2-k\z' — 77 I 2tt\z' — I/ 77 1 


(7) 


which is equal to 


2-7T z" 1 - z’(r\ + i) + 1 


( 8 ) 


Since 0 is placed on the unity circle, the trigonometric identity that describes the unit circle is 0 = 
cos 9 + i sin 9, resulting in 


At the trailing-edge cos 9 


_ r f 

9 2n y £ — cos9 J 

1 , thus the induced velocity at the trailing-edge is 


^&TE 


[ 1 +1 
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(9) 


( 10 ) 


According to the Kutta condition, the total circulation around an airfoil is such that at any instance 
the flow velocity is tangential to the trailing-edge. Therefore, to meet the this condition the tangential 
velocity at the trailing-edge is subtracted from Eq.([9]), thus the total tangential velocity on the airfoil is 


T' ( /7+iY . r' / i - cost) \ / j+I 

27t y£ — cosQ y £ — 1 J 27t\£ — cos9) y £ — 1 


( 11 ) 


The relationship between the velocity vg and the vorticity distribution over the airfoil 7 ( 2 ) according to 
thin airfoil theory is given by the formula y( x ) = —2 vg/ sin0. Further, using the trigonometric relationship 
x = cos 8 and \/l — x 2 = sin 0 , the effect of induced vorticity from a single vortex at a point £ in the wake 
can be written as 


7 i (x,t) = 


/£ + 1 1-x 


( 12 ) 


7 r(£ - x) V £ — 1 V 1 + x 

where ' denote single vorticity element. The effect of the single element of vorticity T can be replaced by 
7 „,(£,t)d£. From Eq.(12) we can derive an expression for the induced vorticity of the entire wake 


. 1 1 — x 

7i(M) = -\/t— 

7T V 1 + X 


iwitt) / £ +1 
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(13) 


Resulting in wake-induced circulation 


r 1 = 
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Unsteady thin airfoil theory 

Karman and Sears [?] applied the principle that, in accordance with the Newton’s second law, the product 
of density and the rate of change of the total momentum at any instance is equal to the lift. In the current 
study, an estimation of the lift due to an unsteady motion is based on the integration of normal pressure 
difference along the chord 


L(x,t) = J Ap(x,f)dx. 


(15) 


The pressure difference A p(x,t) in terms of the chordwise vorticity distribution 7 a (x,f) is expressed by 
the unsteady Bernoulli equation [?] 

d [ x 

Ap{x,t) = pUi a {x,i) + p-^ J 7 Q (x 0 ,t)dx 0 . (16) 

Thus, the lift can be written as 

L(x,t) = pU J j a (x,t) dx + pj IJ 7 a (xo, t) da’o dx. (17) 

Integration by parts of the second integral on the right hand side of Eq. yields 


r 1 d f x 

Pj -Q t J la{xo,t) dx 0 dx = 


dt 


xj la{xo,t) dx 0 


1 

-1 


7 o(x,f)x dx 


where it is recognized that 


X 7a(xo,t)dx 0 


-1 


= J 7 a(M) dx. 


The lift terms in Eq. (17) can then be rearranged as 


f 1 d f 1 d f 1 

L=pU J 70 (x,t)dx-p— J x 7 o(x, t) dx + p— J 7 0 (x,t)dx 

/ 1 d f 1 d f 1 

+ pU J 7 i(x, t) dx — p— J X 7 i(x, t) dx + p— J 7 i(x,t)dx. 


The first term in Eq. (20) 


Lq = pU J 7 o(x, t) dx 


(18) 


(19) 


( 20 ) 


( 21 ) 


is the quasi-steady Kutta-Joukowski lift component. In permanently maintained flow conditions this 
would be the only lift component. In unsteady flow conditions the quasi-steady lift component only 
partially contributes to the total lift and it is determined by evaluating instantaneous angle of attack. 


The second term in Eq.(20) 


d f 1 

L 1 =-p—J X 7 o(x, t) dx 


(22) 
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is the apparent (or added mass) lift component that accounts for the reaction due to the mass of fluid 
directly accelerated by the wing. Following the von-Karman and Sears discussion on time derivative, the 
fifth term of Eq. (20) is equal to 

d 


dt 


/ I O pOO Q p oo 

xo/i(x,t) dx = p— Ji d£ + p—J^ (\/C 2 - 1 -£hw(€,t) d£ 

= 4/ 7 ” ({ ’ ()d f + ov/ (J\ 


ls ±1 


- 1 - 


= P§- t j™i w (U)dt; + P u 7w( ^ } 






(23) 


It should be noted, that the third term of equation Eq. (|23|) is the wake-induced lift 

L 2 =pU 






d£- 


(24) 


Evaluation of L 2 term requires either an assumption about the unsteady motion, keeping track of 
the shed vorticity into the wake, or including the shed vorticity through a convolution integral. The 
forth term in Eq. (201 is cancelled with the second term of equation Eq. (23). Adding the third and sixth 
terms in Eq. (20) with the first term in Eq. (23) results in the Kelvin theorem (see Eq. ®» , therefore, this 
summation is zero. These leads to the time dependent lift (Eq.(20|) that can be written as the sum of 
three terms 


L = pUT 0 -p^J Z 70 (M) dx + puJi d£. 


(25) 


Eq.(25) is the result of von-Karman and Sears [?]. The first term is the quasi-steady lift Lq produced 


by the bound vorticity. The second term L\ represents the apparent (added) mass contribution to lift 
component and it results from the inertia of the fluid moving with the lifting surface. The third term L 2 is 
the induced lift component that produced by the wake vorticity. It should be noted that the contribution 
to the lift due to wake-induced vorticity 71 (x,t) is cancelled out in the derivation of the equations. 

This result coincides with Theodorsen [?], who suggested to divide the time dependent lift into 
circulatory and non-circulatory components, namely Lc and LpjCy respectively. The non-circulatory 
lift can be referred as the virtual mass effect or the acceleration reaction term [?] and the circulatory 
lift component is generated from the vortical flow around the lifting surface. We should point out that 
the non-circulatory lift term, LpiCy which is related to the added mass lift, is identical to the L\ term 


presented by von-Karman and Sears [?] (see Eq.(22)), i.e. L^c = The circulatory lift term Lc is equal 
to the sum of the two remaining lift components, these are the quasi-steady lift and the wake-induced lift, 
i.e. Lc = Lq + L 2 . 


Experimental Apparatus 

Wind tunnel 

The experiments reported herein conducted in a hypobaric climatic wind tunnel at the Advanced Facility 
for Avian Research (AFAR) at the University of Western Ontario. A detailed description of the wind 
tunnel, the experimental technique and the bird is given by Ben-Gida et al. [?] and Kirchhefer et al. [?]. 
Herein we provide a short description, for brevity. The wind tunnel is closed loop type with an octagonal 
test section. The cross-sectional area is 1.2 m 2 , preceded by a 2.5 : 1 contraction ratio. The width, 
height and length of the test section are 1.5 m, lm, and 2 m, respectively. The control of speed, pressure, 
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temperature, and humidity in the wind tunnel enables to simulate flight conditions at high altitudes as 
experienced by birds during long distanced migratory conditions. The bird is introduced into the test 
section through a partition that is located between the downstream end of the test section and the diffuser. 
The turbulence intensity measured by the hot-wires was lower than 0.3% over the entire test section with 
a uniformity of 0.5%. A fine size net was placed at the upstream end of the test section to prevent the bird 
from entering the contraction area. The flight conditions reported in this work correspond to atmospheric 
static pressure, a temperature of 15°C, and relative humidity of 80%. 

The Bird - European Starling 

The wake measurements (as illustrated in Figure [2| were sampled from a European starling (Sturnus 
vulgaris) that was trained to fly in the AFAR wind tunnel. The bird’s wings had an average chord 
of c = 6 cm, a maximum wingspan of b = 38.2 cm ( b sem i = 19.1cm) and an aspect ratio (wingspan 
squared divided by the wings lifting area) of 6.4. The wind speed was set to = 12 m/s. The wingbeat 
frequency, /, was 13.3 Hz on average, and the average peak-to-peak wingtip vertical amplitude, A 1 was 
28 cm. These quantities correspond to a chord-based Reynolds number of 4.8 x 10 4 , a Strouhal number, 
St = Af /Uoo = 0.3, and a reduced frequency, k = nfc/Uoo = 0.2. The birds mass was 78 g and a lateral 
body width of 4 cm. Specially designed safety goggles (Yamamoto Cogaku Co. model YL600) were 
adjusted to the bird while flying at the tunnel. In addition, a collection of optoisolators operated by six 
infrared transceivers were integrated into the PIV system in order to prevent direct contact between the 
bird and the laser sheet. The optoisolators triggered the laser only when the bird was flying upstream 
further from the PIV field of view. All animal care and procedures were approved by the University of 
Western Ontario Animal Use Sub-Committee (protocols 2006-011, 2010-216). 

Long duration time resolved PIV 

Flow measurements were taken using the long-duration time-resolved PIV system developed by Taylor 
et al. [?]. The PIV system consists of a 80 W double-head diode-pumped Q-switched Nd:YLF laser at a 
wavelength of 527 nm and two CMOS cameras (Photron FASTCAM-1024PCI) with spatial resolution of 
1024 x 1024pixel 2 at a sampling rate of 1000 Hz. The PIV system is capable of acquiring image pairs at 
500 Hz using two cameras for a continuous period of 20 minutes. Olive oil aerosol particles, 1 /im in size 
on average [?] were introduced into the wind tunnel using a Laskin nozzle from the downstream end of 
the test section so that it did not cause a disturbance to the flow in the test section or to the bird. The 
system is designed to work either in a 2D or Stereo mode. Herein, we used one camera for the PIV whilst 
the other camera was used for measuring the wingbeat kinematics simultaneously with the PIV. The 
PIV camera’s field of view was 12 x 12cm 2 corresponding to 2c by 2c. The velocity fields were computed 
using OpenPIV [?] using 32pixel 2 interrogation windows with 50% overlap, yielding a spatial resolution 
of 32 vectors per average chord, equal to 1.8 vectors per millimetre. In the current experiments, 4,600 
vector maps were recorded, and out of this dataset 650 vector maps contained features of the near wake 
behind the starling’s wing. The measured wake was located 4 wing chord lengths behind the right wing. 
The wake was sampled in the streamwise-normal plane at 2ms intervals (500Hz), therefore, both the 
downstroke and the upstroke phases were temporally resolved. 

Kinematic measurements 

To relate the wake measurements to the kinematic motion of the bird’s wings, an analysis of the kinematic 
motion has been undertaken. The simultaneous measurements also enables a point of comparison between 
the estimation of lift through the unsteady thin airfoil theory and lift calculation using the momentum 
equation for viscous flows. The field of view by the CMOS camera is 9c by 9c (corresponding to 54 x 
54cm 2 ). Figure [ 2 ] depicts a sample image of the Starling flying in the wind-tunnel as captured by the 


camera. The box marked with the ‘PIV’ label indicates the location of the measured velocity fields from 
the PIV system. In addition, a floor-mounted camera operating at 60 Hz was used to record the spanwise 
position of the bird with respect to the laser sheet illumination. Therefore, these images provided the 
identification of the measured PIV plane and its location in respect to the wing; so that, the wake velocity 
field associated with the spanwise location across the wing is 0.15&/2 from the root. The floor-mounted 
camera was not synchronized with the PIV or the kinematic measurements, so the two time histories 
were synchronized manually based on the presence of the laser light in the images. Once synchronized, 
spanwise positions were assigned to the wake data acquired at 500 Hz based on interpolation from the 
simultaneously recorded spanwise positions. 

Error estimation 

An error analysis based on the root sum of squares method has been applied to the velocity data and 
the wing kinematics. The errors were estimated as: 2.5% for the instantaneous velocity values, 12% for 
the instantaneous vorticity and 3% for the lift values [?]. The error introduced in the kinematic analysis 
resulted from the spatial resolution of the image and the lens distortion leading to an estimated error of 
5% in the wing displacements. 


Results and Discussions 

The PIV flow field measurements and the bird’s kinematics were each analysed separately in order to 
estimate the time dependent component of the lift generated by the flapping motion of the starling using 
the inviscid and the viscid approaches. Linear lift theory (see the theoretical modelling section) was used 
to estimate the lift from the bird’s kinematics, whereas a viscous flow theory, derived by Wu [?] and 
applied by Panda [?], was implemented to estimate the lift from the near wake flow fields measured by 
the PIV. The present work considers a comparison between the linear theories [?,?] and the viscous flow 
theory [?], for the lift generated by the starling. Both approaches will emphasize the role of the time 
dependent lift components. 

The bird kinematics and the near wake velocity field were captured simultaneously while the starling 
was flying in the tunnel and flapping its wings continuously. During flapping, birds generate lift and 
thrust. The lift is necessary to support the bird’s body weight and the thrust is required to overcome 
the drag. The data presented herein correspond to no-maneuver and no-acceleration conditions. The 
spatial location of the bird’s body at the beginning of the downstroke and upstroke phases of flight are 
shown in Figure [3] for three consecutive wingbeats. It can be seen that the bird does not accelerate in the 
streamwise or vertical directions [?]. Hence, the following kinematic analysis can be performed assuming 
negligible acceleration. 

Estimation of time-dependent lift from the bird kinematics 

As presented in the theoretical modelling section, the linear theories were derived within the framework of 
potential theory, which assumes inviscid fluid with small disturbances and a plane wake. One can use 
such theories to estimate the time dependent lift components from the kinematics of a wing section with 
relatively good precision [?]. We choose to use the guidelines provided by Leishman [?] and estimate the 
time dependent lift from the bird’s kinematics. 

Using the aforesaid unsteady thin airfoil theory, we estimate the lift generated by the flapping wings 
motion of the European starling, as captured through the wings kinematic images. For simplicity, we 
describe the kinematics of a flapping wing by a pure two-dimensional plunging motion, which involves a 
heaving up and down of the wing section that results in a variation of the effective angle of attack [?]. In 
such motion the variation of the vertical displacement with time can be described as follows 
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hf(t) = hf 0 cos (ut) (26) 

where hf 0 is the plunging motion amplitude and u> is the angular velocity. For a flapping wing the plunging 
amplitude hf 0 is a function of the distance from the shoulder joint. Assuming the elastic deformation 
of the wing is negligible we can describe the plunging amplitude as linearly increasing function towards 
the wing tip. Thus, the vertical displacement amplitude at spanwise distance rj = 0.156/2 from the wing 
root can be reduced to hf 0 = ?7cos(0), where <j> is the wing tip angle, as depicted in Figure[4j Kinematic 
images of the starling [?] depicted the wings as they oscillate in a periodic manner, where the range of the 
angular positions is —55° < <fi < 19° 

Due to the unique bone and muscle structure of the bird’s wing, during flapping flight the inner 
part of the wing experiences less twisting motion than the outer part which accounts for most of the 
thrust production [?]. Hence, the variation of the local angles of attack at the inner part of the wing is 
small compared to those at the outer part. Therefore, we can describe the effective angle of attack as a 
result of the horizontal free-stream velocity and the vertical velocity component due to plunging motion, 
a e = tan ~ 1 (h(t)/U oa ). By assuming small angles of attack we can simplify the effective angle of attack 
to a e ~ (h(t)/Uoo). Consequently, the quasi-steady lift component can be estimated by the thin airfoil 
theory [?], where the non-dimensional lift coefficient is a function of effective angle of attack and the 
corresponding lift component is equal to 


L 0 = npU 2 c 


h(t) 

U 


(27) 


Following the unsteady thin airfoil theory [?,?], the added-mass lift component is a result of flow 
acceleration, and thus arises from the unsteady term in the Bernoulli equation that accounts for the 
pressure force required to accelerate the fluid in the the vicinity of the wing. For the wing section moving 
normal to its surface at velocity v(t), the non-circulatory fluid force acting on the surface is equal to the 
product of apparent mass and acceleration. Thus, a body moving in an unsteady motion must overcome 
acceleration in addition to its own inertial force. Therefore, the apparent mass (or non-circulatory [?]) lift 
component can be estimated from the kinematic motion accordingly 


L\ 


= npU 2 


4 


Hi) 

U' 2 


(28) 


where h(t) is the time derivative of vertical displacement [?]. 

Here, the effect of the wake-induced lift component L 2 is determined by assuming harmonic motion [?] 
at a frequency w, yielding j w (£,t) = ge lUJ< - t ^^^ u K Using Theodorsen’s function C(k), which accounts for 
the effect of the shed vortices on the unsteady aerodynamic loads, we can calculate the wake-induced lift 
component as follows 

L 2 (t ) = (C(k) - 1)L 0 (29) 

where the definition of Theodorsen’s function is [?] 


C(k) 


H[ 2 \k) 

H[ 2) {k)+iH^\k)' 


(30) 


( 2 ) 

Here H n = J v — iY v is the Hankel function of the reduced frequency k, where J v and Y v are Bessel 
functions of the first and second kind, respectively. The lift reduction function C(k) falls gradually to 
value of 0.5 as k goes to infinity. The effect of C(k) as producer of phase lag takes over very quickly, 
where the maximum rotation of the vector occurs around k = 0.2. It should be noted that the reduced 
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frequency values used by many small passerines, such as the European starling, lays well in the range of 
0.1 < k < 0.3 [?]. Apparently, these small passerines fly at the region where the effect of the lift reduction 
function is the strongest. 

Figure [5] shows the time variation of the three lift components (Lo, Ti and L 2 ). The sum of the 
three components correspond to the total time dependent lift L generated during the flapping motion. 
It appears that the non-circulatory (or added mass) contribution to lift is the smallest among the three 
components. Yet, this contribution is not negligible and in fact is equal to about half of the the lift 
generated by the induced-wake. Both components are significantly smaller compared to the quasi-steady 
lift component. According to the harmonic assumption the wake-induced lift is in anti-phase to the 
quasi-steady lift component. The generation of the circulatory lift components comprises of two terms 
that in fact are counter to each other. The theory explains the negative work done by the induced vorticity 
during the upstroke phase as the wing approaches maximum value of the lift. Rival et al. [?] utilized 
similar principles to study the effect of leading-edge vortex on the formation of dynamic stall over an 
airfoil. 

Estimation of the circulatory time-dependent lift component from the bird’s 
near wake vorticity field 

The unsteady thin airfoil theory is a useful tool that provides a good approximation of the time de¬ 
pendent aerodynamic loads. Nevertheless, the theory, which bounds to two-dimensional inviscid flows, 
underestimates loads in complex flows where viscosity cannot be ignored. Therefore, another approach is 
needed. One of the methods for estimating aerodynamic loads is based on the flow field at the wake region. 
PIV provides high resolution spatial data with sufficient accuracy that enables the estimation of such 
loads from the wake of bluff bodies [?]. However, despite many advances in the current state-of-the-art 
in experimental diagnostics, practical application of PIV to estimate time dependent forces from wake 
flow-field measurements are challenging. These efforts are limited to 2D planes. This limitation is mainly 
due to the fact that the 3D flow-field measurements are restricted by relatively small volume size and low 
Reynolds number flows. Thus, the most practical approach to estimate the time dependent lift component 
is from 2D plane measurements. 

Within the avian research community, the most common 2D approach is concerned with measurements 
of the vorticity field in the far wake Trefftz plane and application of the Kutta-Joukowski quasi-steady 
theorem in order to estimate the lift. This approach is based on the classical assumption that the vortex 
lines behind a lifting surface roll-up when they propagate downstream into the wake, and they bundle 
into tip vortices. Thus, the far wake is dominated by the tip vortices. This approach is appealing as the 
entire wake structures can be captured by a single plane, provided these measurements are acquired in 
the far wake. However, measurements conducted in the Trefftz plane are highly inaccurate due to plane 
normal velocity component and may lead to significant errors that are hard to ignore [?]. Furthermore, 
these measurements allow estimation of only the mean quasi-steady total lift and not the time dependent 
evolution of the lift. 

The second 2D approach is concerned with measurements of the flow field in the streamwise plane at 
the near or far wake. A brief summery of the PIV measurement acquired in the wake of freely Hying birds 
can be found in figure 1 at Kirchhefer et al. 2013 [?]. Amongst which are flow measurements in the wake 
of Thrush Nightingale [?], Robin [?], Swift [?] and bats [?,?]). PIV measurements in the streamwise plane 
are considered to be accurate, with some errors related to the spanwise velocity component [?]. As it 
has been indicated previously, one of the first applications of PIV technique to estimate the lift from the 
wake of freely flying bird is attributed to Spedding et. al. [?]. In this work, the lift was estimated based 
on quasi-steady Kutta-Joukowski thin airfoil theory. This simplified approach, which has been followed 
by many other researchers, neglects the effects of added mass and wake-induced vorticity on the time 
dependent lift components. 
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In the current work, the near wake flow-fields were captured simultaneously to the bird’s kinematic 
motion, shown in Figure [3] and discussed in the previous section, thus, allowing one to relate the wake 
flow-field structures to the bird’s kinematic motion. The PIV measurements were taken at the inner part 
of the right wing (from the bird’s perspective), at a spanwise distance of z = 0.15&/2 from the wing root. 
In order to shed light on the wake structures that manifest the bird’s lift, a visualization of the entire 
wingbeat during a single flapping cycle is performed by generating a wake composite image from multiple 
PIV realizations. A similar approach was first applied by Spedding et al. [?] in which PIV measurements 
(from separate wingbeat cycles) were arranged to represent a complete and representative wavelengths of 
the wake. 

The wake composite is formed by plotting sequential PIV realizations, each image is offset to one- 
another in the streamwise direction. The offset of the PIV images is calculated as Ua^At-n. The generation 
of a wake composite provides a useful tool for observing the time-series of measurements representing 
the wake of a wingbeat cycle. The procedure was performed using the PIV flow-fields collected at a 
sampling rate of 500 Hz that is significantly higher than the bird’s flapping frequency of 13 Hz. Therefore, 
a pattern of vorticity appearing in one frame also appears in the consecutive frame - only phase-shifted. 
The wake structures that appear ‘downstream’ in the wake composite image happen earlier in time, while 
the structures that appear ‘upstream’ in the composite actually happen later in time. In a sense, the 
generation of the wake composite image invokes Taylors hypothesis [?] in which the characteristics of the 
flow are advected through the field of view, where the offset of one image to the next is based on the free 
stream speed. It should be noted that the typical offset of U^At ■ n between images is ~ 0.4c and an 
instantaneous PIV measurement has a spatial dimension of 2c. Therefore, at any location in the wake 
composite image, there are several overlapping images that can be used to ascertain the instantaneous 
wake characteristics over the streamwise distance of 2c to compare with the wake composite at the same 
location. 

The wake features are shown through fluctuating velocity and vorticity fields, where the spanwise 
vorticity is defined as follows 


Uz{t) = 


dv 

dx 


du 

dy 


(31) 


and is evaluated directly from the PIV flow fields using a least squares differentiation scheme. Here u and 
v are streamwise and transverse velocity components, respectively. 

Estimation of time dependent lift generated during the flapping motion of the starling is evaluated from 
the near wake velocity maps by utilizing the viscid approach derived by Wu [?] based on the Navier-Stokes 
equation [?,?]. The generalized formulation that conveniently describes the aerodynamic forces exerted 
by a fluid on a solid body immersed in, and moving relative to the fluid, is equal to inertial force due to 
the mass displaced by the solid body and a term proportional to the time of change of the first moment of 
the vorticity field [?,?], as follows: 


L{t) = ~ p dt 


xui z (t) dxdy 


m 


,dU 

dt 


(32) 


where m' is the mass of the fluid displaced by the solid body. One can immediately recognize that the 
second term in Eq.(32) is the added mass lift component, which correpsonds to the L\ lift component 


in the von-Karman and Sears notation. The term ff xco z dxdy represents the first r-moment of the 
vorticity. The equation derived by Wu [?] is based on the principle that if the vorticity distribution over 
the entire flow field were known the force could be evaluated accurately. Although the lift terms may, for 
utility and convenience, be divided farther into Lq and L 2 components, such division is to some extent 
arbitrary. According to this approach, the circulatory lift L c (t) is equal to the time rate of change of the 
first moment of the vorticity field. By applying the Taylor hypothesis, dx = U c dt, one can transform 
the spatial derivative into a temporal one. In the unsteady thin airfoil terminology the circulatory lift 
component is only a portion of the total lift that acting on the flapping wing. Since at the beginning 
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of the flapping cycle the lift is unknown we refer to the estimated lift component as an increment in 
the circulatory lift that is generated from the beginning of the cycle, thus equal to A L c (t) and can be 
expressed as 


A L c (t) = pUoc J U c ((t)dt 


(33) 


In order to estimate the circulatory lift AL c (t) from Eq.(33) one needs to acquire information regarding 
the vorticity flux £(t) in the near wake 


C(t) = J u z {t)dy 


(34) 


The vorticity flux, defined by Eq.(34), is estimated for each individual vector map as function of time. 


The calculated vorticity flux corresponds to the spanwise vorticity component and is integrated over a 
selected region in each vector map. The selected region covers the wake features that are observed in 
figures 6(b) 7(b)| and |8(b)~[ Figure 6(a) demonstrates the changes in the circulatory lift component 
as it evolves over time, calculated based on Eq. (33) over a single wingbeat cycle. This calculation was 


performed for three different wingbeat cycles. The curve represents the cumulative lift over one wingbeat 
cycle, starting from right to left as the bird is moving from right to left in respect to the coordinate 
system. Overall, the lift accumulates positively during the downstroke phase whilst negative accumulation 
is depicted during the upstroke phase. One can deduce that the circulatory lift has a positive net effect 
during the downstroke phase, which is in agreement with former work [?,?]. During the upstroke phase, 
it appears that the circulatory lift is decreasing; this implies that during this phase, the bird is losing 
energy through the unsteady mechanism. Furthermore, the presence of cumulative negative circulatory 
lift during the upstroke phase marks the energy that the bird has to invest in order to bring the wing 
back to the downstroke phase to generate lift, again. 

In the case of flapping flight of natural free-flyers, the wings’ motion is extremely complicated and 
it comprises from a complex flapping motion (changes in effective angle of attack), wing deformation 
(the bird stretches or bending its wing) and substantial three dimensional motion [?], thus resulting in a 
complex wake vorticity system [?]. As was shown earlier by the unsteady thin airfoil theory the vorticity 
shed into the wake continuously and affect the total circulation around a lifting surface [?,?]. Therefore, 
as demonstrated in figures [6j [TJ and [8j the evolution of the lift over a single wingbeat cycle should be 
considered even for the case of power estimates where it is shown that while on average the lift should 
be equal to the bird’s weight, the time dependent variations of the lift from this value might provide a 
plausible argument to a more efficient flight. 

Pennyquick [?] suggested to use the quasi-steady approach when applied to estimating lift in bird 
flight. In his work, power was estimated based on kinematic analysis of flying birds and some assumptions 
related to drag and lift. The lift was assumed to be equal to weight, as any body that is aloft and in 
equilibrium. Following this approach, Tucker [?] revisited this argument and refined it to consider the 
flapping motion of the wings. Based on these works, Rayner [?] proposed a mathematical model for 
lift estimation from the wake of a flying bird in various flight modes. This model is well supported in 
the literature and it is conceptually accepted that unsteady mechanisms are of minor importance since 
variation of wing pitch and circulation are not producing thrust [?]. Our results, on the other hand, 
demonstrate that the unsteady mechanisms indeed play a significant role in the generation of lift. Whilst, 
thrust is not generated by lift, it is the energy that is required to keep the bird aloft that is impacted by 
the lift mechanism, e.g.: with less power required to generate lift, more power can be directed towards 
propulsion. 

The argument that during the steady phase of flight the bird’s weight must be balanced by an equal 
amount of lift [?] is obviously valid for an average lift generated during a wingbeat cycle. In the current 
study the bird’s weight was equal to 2 N/m. Following this argument, for steady flapping flight, the 
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starling is required to generate an equal amount of lift force. The three cycles show similar trends, where 
the lift variation over each wingbeat cycle is about ±2 N/m. As mentioned above, these lift values are 
the time variation of the circulatory lift component. Thus, in order to obtain the actual total lift force 
produced by the starling, one needs to add these fluctuations to the bird’s weight that represents the 
averaged lift over a wingbeat cycle. By adding the lift fluctuations with the starling’s weight we find that 
during the flapping cycle the starling can produce up to 4 N/m of lift force, which is also equal to twice 
of its body weight. These maximum lift values are generated mostly during the downstroke to upstroke 
transition phase, as depicted in figure 6(b)[ 7(b) and |8(b)| which is in good agreement with Kirchhcfer et 
al. 2013 [?] who showed particular flow structures (termed: ’double branch features’) occurring during 
the downstroke to upstroke transition. In addition, one can observe that the the minimum lift values 
are produced during the transition from upstroke to downstroke phase. The physical argument to this 
observation is that during the upstroke to downstroke transition the starling folds its wings, causing them 
to stop acting as lifting surfaces, and thus generates almost no lift. The aforesaid results imply that the 
usage of quasi-steady lift theory might be underestimating the lift that a bird is actually generates during 
flapping flight. 


Conclusions 

The objective of this study was to evaluate the effect of flapping motion on the generation of unsteady 
lift. Therefore, freely flying European starling’s kinematics and near wake flow fields were acquired 
simultaneously at the AFAR facility. The bird’s kinematics were measured with high speed imaging 
system, whilst the near wake was acquired with long duration time resolved PIV. 

To estimate the time dependent lift we have applied unsteady thin airfoil theory as developed by 
Theodorsen [?] and von-Karman and Sears [?]. The theory addresses the various mechanisms that 
contribute to the time dependent lift components: the quasi-steady, added mass and wake-induced 
vorticity. Using these terms, we have estimated the lift from the wingbeat kinematics and demonstrated 
the contribution of each one of the terms to the total lift. 

The theory assumes a planar wake and a trailing-edge Kutta condition. Thus it excludes wake roll-up, 
convection of large separations over the wing section, boundary-layer separation, large laminar separation 
bubbles, leading-edge and trailing-edge vortices, three-dimensional effects and so forth. However, these 
effects are significant at low Reynolds numbers and thus lift estimated with the unsteady thin airfoil 
theory only partial explains the complex flow physics. 

The aerodynamic forces estimated from the wake flow field measurements as acquired with PIV provide 
a more accurate estimation of the aerodynamic forces. The equations that provide the basis to determine 
the aerodynamic forces were derived by Wu [?]. As discussed by Panda [?] only the circulatory lift 
component can be estimated from the wake measurements. In this work, three wingbeat cycles were 
analysed. The near wake behind the bird’s wing was spatially reconstructed from the temporal data. The 
reconstructed wake field is spanning over 20 chord lengths downstream and depicts the flow features in 
the wake. The circulatory lift component over the wingbeat cycle follows the flow features as shown by 
the reconstructed wake. 

The evolution of the circulatory lift presents a negative contribution during the upstroke phase, whilst 
during the downstroke phase a positive contribution is observed. In addition, we observed that the 
downstroke phase is shorter, compared to the upstroke phase. This asymmetrical pattern is essential 
for the production of the high lift impulse that supports the bird’s weight. We conclude that the bird 
prefer to generate high lift values with a short downstroke phase than moderate lift values with a longer 
downstroke phase. 

The variation of the lift over the wingbeat cycle emphasizes its contribution to the total lift and its 
role in power estimations. We suggest that the time dependant circulatory lift component cannot be 
assumed negligible and should be considered when estimating lift or power of birds in flapping motion. 
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Figure 1. Auxiliary diagram showing notation employed. 
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0.6 m 


Figure 2. The large image shows the kinematic camera field of view and the small window 
marked ‘PIV’ is the PIV camera field of view. 
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(c) wingbeat 2-upstroke (d) wingbeat 2—downstroke 



(e) wingbeat 3-upstroke (f) wingbeat 3-downstroke 


Figure 3. Side view of the European Starling at the wind tunnel. The white line is placed to 
provide spatial reference of the bird’s body. The white line is inclined at 8.8° with the free-stream 
velocity. The left and right images correspond to beginning of the downstroke and upstroke, respectively. 
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Figure 4. Auxiliary diagram of the flapping motion with the definition of the flapping 
angle. View in the direction of the flow. 
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Figure 5. Time dependant components of unsteady lift terms, as estimated from the wing 
kinematics. The terms Lq, Li, L 2 , are estimated from Eqs. 27j 28 29j respectively. The term 
L = Lq + L\ + L 2 is the total lift. 
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Figure 6. Estimation of circulatory lift component based on wingbeat number 1. 


(a) The circulatory lift was estimation based on eq.(33). The grey area indicates downstroke flapping 
phase, (b) Reconstruction of the starling’s wake vorticity as thought the bird flies from right to left. 
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Figure 7. Estimation of circulatory lift component based on wingbeat number 2. 


(a) The circulatory lift was estimation based on eq.(33). The grey area indicates downstroke flapping 


phase, (b) Reconstruction of the starling’s wake vorticity as thought the bird flies from right to left. 
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Figure 8. Estimation of circulatory lift component based on wingbeat number 3. 

| (a) | The circulatory lift was estimation based on eq. (33). The grey area indicates downstroke flapping 
phase, (b) Reconstruction of the starling’s wake vorticity as thought the bird flies from right to left. 














